#include "../includes/need-a-name.hh"
#include "includes/initcollecter-sound-wave.hh" 
#include <string>

InitCollecterSoundWave::InitCollecterSoundWave(){
  set_parameter("wid1",0.0);
  set_parameter("wid2",0.0);
  set_parameter("wid3",0.0);
  set_parameter("r",0.0);
  set_parameter("k1",0.0);
  set_parameter("k2",0.0);
  set_parameter("k3",0.0);
  set_parameter("cen1",0.0);
  set_parameter("cen2",0.0);
  set_parameter("cen3",0.0);
  set_parameter("rho_rectangle",0.0);
  set_parameter("rho_circle",0.0);
  set_parameter("rho_plane",0.0);    
}
InitCollecterSoundWave::~InitCollecterSoundWave(){
}

void InitCollecterSoundWave::initial_function(double* location, double pillarlen, int currentdim, int index){
  double componenttmp;  
  double basisInside[memory::dgorder*memory::intergnumber];
  double locationtmp[gridsetting::dimension];  
  if(memory::dgflag==0){
    for(int tmpi=0; tmpi<memory::phynumber;tmpi++){
      if(para["rho_rectangle"]==1.0)
	memory::variable[tmpi][index]=rho_rectangle(location,tmpi);
      if(para["rho_circle"]==1.0)
	memory::variable[tmpi][index]=rho_circle(location,tmpi);
      if(para["rho_plane"]==1.0)
	memory::variable[tmpi][index]=rho_plane(location,tmpi);
    }
  }
  if(memory::dgflag==1){
    basis.getBasis(basisInside);
    for(int tmpl=0; tmpl<gridsetting::dimension; tmpl++)
      locationtmp[tmpl]=location[tmpl];
    for(int tmpi=0; tmpi<memory::phynumber;tmpi++)
      for(int tmpj=0; tmpj<memory::dgorder;tmpj++){
    	componenttmp=0.0;
    	for(int tmpk=0; tmpk<memory::intergnumber; tmpk++){
    	  locationtmp[currentdim]=location[currentdim]+(-0.5+tmpk/(memory::intergnumber-1.0))*pillarlen;
	  if(para["rho_rectangle"]==1.0)
	    componenttmp+=(2.0/memory::intergnumber)*(2.0*tmpj+1.0)/2.0*rho_rectangle(locationtmp,tmpi)*basisInside[tmpj*memory::intergnumber+tmpk];
	  if(para["rho_circle"]==1.0)
	    componenttmp+=(2.0/memory::intergnumber)*(2.0*tmpj+1.0)/2.0*rho_circle(locationtmp,tmpi)*basisInside[tmpj*memory::intergnumber+tmpk];
	  if(para["rho_plane"]==1.0)
	    componenttmp+=(2.0/memory::intergnumber)*(2.0*tmpj+1.0)/2.0*rho_plane(locationtmp,tmpi)*basisInside[tmpj*memory::intergnumber+tmpk];
	}
	memory::variable[tmpj*memory::phynumber+tmpi][index]=componenttmp;
      }
    for(int tmpi=0; tmpi<memory::phynumber;tmpi++){
      if(para["rho_rectangle"]==1.0)
    	memory::variable[memory::dgorder*memory::phynumber+tmpi][index]=rho_rectangle(location,tmpi);
      if(para["rho_circle"]==1.0)
    	memory::variable[memory::dgorder*memory::phynumber+tmpi][index]=rho_circle(location,tmpi);
      if(para["rho_plane"]==1.0)
    	memory::variable[memory::dgorder*memory::phynumber+tmpi][index]=rho_plane(location,tmpi);
    }
  }
}

double InitCollecterSoundWave::rho_rectangle(double* location, int phyi){
  double cen[gridsetting::dimension];
  double wid[gridsetting::dimension];
  double tmp=0.0;
  if(gridsetting::dimension>=1){
    cen[0]=para["cen1"];
    wid[0]=para["wid1"];
  }
  if(gridsetting::dimension>=2){
    cen[1]=para["cen2"];
    wid[1]=para["wid2"];
  }
  if(gridsetting::dimension==3){
    cen[2]=para["cen3"];
    wid[2]=para["wid3"];
  }
  if(phyi==0){
    tmp=1.0;
    for(int tmpi=0; tmpi<gridsetting::dimension;tmpi++)
      tmp=tmp*(std::abs(location[tmpi]-cen[tmpi])<wid[tmpi]? 1.0 : 0.0);  
  }
  if(phyi>=1){
    tmp=0.0;
  }
  return tmp;
};

double InitCollecterSoundWave::rho_circle(double* location, int phyi){
  double cen[gridsetting::dimension];
  double radius;
  double distance=0;
  double tmp=0.0;
  if(gridsetting::dimension>=1)
    cen[0]=para["cen1"];
  if(gridsetting::dimension>=2)
    cen[1]=para["cen2"];
  if(gridsetting::dimension==3)
    cen[2]=para["cen3"];
  radius=para["r"];
  for(int tmpi=0; tmpi<gridsetting::dimension;tmpi++)
    distance=distance+std::pow(location[tmpi]-cen[tmpi],2);
  distance=std::sqrt(distance);

  if(phyi==0)
    tmp=(distance<radius? 1.0 : 0.0);
  if(phyi>=1)
    tmp=0.0;
  return tmp;
};

double InitCollecterSoundWave::rho_plane(double* location, int phyi){
  double cen[gridsetting::dimension];
  double k[gridsetting::dimension];
  double tmp=0.0;
  if(gridsetting::dimension>=1){
    cen[0]=para["cen1"];
    k[0]=para["k1"];
  }
  if(gridsetting::dimension>=2){
    cen[1]=para["cen2"];
    k[1]=para["k2"];
  }
  if(gridsetting::dimension==3){
    cen[2]=para["cen3"];
    k[2]=para["k3"];    
  }
  double dotwithnormal=0;
  for(int tmpi=0; tmpi<gridsetting::dimension;tmpi++)
    dotwithnormal=dotwithnormal+k[tmpi]*(location[tmpi]-cen[tmpi]);

  if(phyi==0)
    tmp=(dotwithnormal>0? 1.0 : 0.1);    
  if(phyi>=1)
    tmp=0.0;
  return tmp;
}
